package br.unb.statistic;

import br.unb.utils.Resource;
import java.io.IOException;

public class SpectrumProcess {

    private double hurst;
    private double p;
    private double q;

    public SpectrumProcess(double hurst) {
        super();
        setHurst(hurst);
    }

    public double getHurst() {
        return hurst;
    }

    public void setHurst(double hurst) {
        double F, G, pi2, pi3;

        this.hurst = hurst;

        pi2 = Math.pow(Math.PI, 2);
        pi3 = Math.pow(Math.PI, 3);
        F = F(this.hurst);
        G = G(this.hurst);

        this.p = (-6 / pi2) * F + (12 / pi3) * G;
        this.q = (4 / Math.PI) * F - (6 / pi2) * G;

    }

    private double F(double H) {
        double ret = 0;
        double a, b;

        for (int i = 3; i < 204; i++) {
            a = 2 * i * Math.PI - Math.PI;
            b = 2 * i * Math.PI + Math.PI;
            ret += (Math.pow(a, -(2 * H)) - Math.pow(b, -(2 * H)));
        }
        ret = ret / (2 * H);
        return ret;
    }

    private double G(double H) {
        double ret = 0;
        double a, b;

        b = -2 * H;
        if (H != 0.5) {
            for (int i = 3; i < 204; i++) {
                a = 2 * i * Math.PI;
                ret += (a * Math.pow((a + Math.PI), b) + a * Math.pow((a - Math.PI), b) - 2 * Math.pow(a, (b + 1))) / (2 * H);
                ret -= (Math.pow((a + Math.PI), b + 1) + Math.pow((a - Math.PI), b + 1) - 2 * Math.pow(a, (b + 1))) / (2 * H - 1);
            }
            return ret;
        } else {
            for (int i = 3; i < 204; i++) {
                a = 2 * i * Math.PI;
                ret += (a / (a + Math.PI) + a / (a - Math.PI) + Math.log(a + Math.PI) + Math.log(a - Math.PI) - 2 - 2 * Math.log(a));
            }
            return ret;
        }
    }

    private double B(double lambda) {
        double a, b, ret = 0;
        double B3 = this.p * lambda + this.q;

        b = (-2 * this.hurst) - 1;
        for (int i = 1; i < 3; i++) {
            a = 2 * i * Math.PI;
            ret += Math.pow((a + lambda), b) + Math.pow((a - lambda), b);
        }
        ret += B3;
        return ret;
    }

    public static double gammaln(double xx) {
        double cof[] = {
            76.180091730000001D, -86.505320330000004D, 24.014098220000001D, -1.231739516D, 0.00120858003D, 5.3638199999999999E-006D
        };
        double x = xx - 1.0D;
        double tmp = x + 5.5D;
        tmp -= (x + 0.5D) * Math.log(tmp);
        double ser = 1.0D;
        for (int j = 0; j <= 5; j++) {
            x++;
            ser += cof[j] / x;
        }

        return -tmp + Math.log(2.5066282746500002D * ser);
    }

    public static double gamma(double a) {
        return Math.exp(gammaln(a));
    }

    public double[] fgnSpectrum(int n) {
        double lambda, a, b, c;
        double A = 2 * Math.sin(Math.PI * this.hurst) * gamma(2 * this.hurst + 1);
        double PowSpectrum[] = new double[n / 2 + 1];
        PowSpectrum[0] = 0;
        for (int i = 1; i <= n / 2; i++) {
            lambda = (2 * i * Math.PI) / n;
            a = A * (1 - Math.cos(lambda));
            b = Math.pow(Math.abs(lambda), (-2 * this.hurst - 1));
            c = B(lambda);
            PowSpectrum[i] = a * (b + c);
        }

        return PowSpectrum;
    }

    //começo
    public double[] fgnSpectrum2(int i) {
        double theta;
        double d1 = this.hurst;
        double d2 = -2D * d1 - 1.0D;
        char c = '\310';
        int j = (int) Math.floor((double) (i) / 2D) + 1;
        double ad[] = new double[j];
        double ad1[] = new double[j];
        double ad2[] = new double[j];
        double d3 = 0.31830988618379069D * gamma(-d2) * Math.sin(3.1415926535897931D * d1);
        double d10 = 6.2831853071795862D / (double) i;
        for (int k = 0; k < j; k++) {
            ad[k] = d10 * (double) (k + 1);
        }

        for (int j1 = 0; j1 < j; j1++) {
            double d5 = 0.0D;
            for (int k1 = 1; k1 <= c; k1++) {
                double d8 = Math.abs(6.2831853071795862D * (double) k1 + ad[j1]);
                double d9 = Math.abs(-6.2831853071795862D * (double) k1 + ad[j1]);
                d5 = d5 + Math.pow(d8, d2) + Math.pow(d9, d2);
            }

            double d7 = Math.pow(ad[j1], d2);
            double d6 = d5 + d7;
            ad2[j1] = (1.0D - Math.cos(ad[j1])) * d3 * d6;
            ad1[j1] = Math.log(ad2[j1]);
        }

        double d4 = 0.0D;
        for (int l = 0; l < j; l++) {
            d4 += ad1[l];
        }

        theta = Math.exp((d4 * 2D) / (double) i);
        for (int i1 = 0; i1 < j; i1++) {
            ad2[i1] = ad2[i1] / theta;
        }

        return ad2;
    }
    //fim

    public static double[] periodogram(double timeDomainRealSeries[], int samples) throws IOException {
        double[] input = new double[2 * samples];
        for (int i = 0; i < samples; i++) {
            input[2 * i] = timeDomainRealSeries[i];
            input[2 * i + 1] = 0;
        }
        double[] fft = SpectrumProcess.fft(input, samples, 1);
        double[] periodogram = new double[fft.length / 4];
        for (int i = 0; i < fft.length / 4; i++) {
            periodogram[i] = ((fft[2 * i] * fft[2 * i]) + (fft[2 * i + 1] * fft[2 * i + 1])) / (2 * Math.PI * samples);
        }
        return periodogram;
    }

    public static double[] fft(double input[], int samples, int isign) throws IOException {

        int i;

        double tempr, tempi, theta, wtemp, wpr, wpi, wr, wi;

        //Verifies if the number of samples is a power of 2
        for (i = 2; i <= 32 && Math.pow(2, i) != (double) samples; i++);

        if (i == 33) {
            throw new IOException(Resource.getString("NUMBER_MUST_BE_A_POWER_OF_2"));
        }

        int n = samples << 1;

        double output[] = new double[n + 1];

        for (i = n; i > 0; i--) {
            output[i] = input[i - 1];
        }

        int j = 1;
        for (i = 1; i < n; i += 2) {
            if (j > i) {
                tempr = output[j];
                output[j] = output[i];
                output[i] = tempr;
                tempr = output[j + 1];
                output[j + 1] = output[i + 1];
                output[i + 1] = tempr;
            }
            int m = n >> 1;
            while (m >= 2 && j > m) {
                j -= m;
                m >>= 1;
            }

            j += m;
        }

        int istep;
        int mmax = 2;
        while (n > mmax) {
            istep = 2 * mmax;
            theta = (2 * Math.PI) / (double) (isign * mmax);
            wtemp = Math.sin(0.5D * theta);
            wpr = -2D * wtemp * wtemp;
            wpi = Math.sin(theta);
            wr = 1.0D;
            wi = 0.0D;
            for (int m = 1; m < mmax; m += 2) {
                for (i = m; i <= n; i += istep) {
                    j = i + mmax;
                    tempr = wr * output[j] - wi * output[j + 1];
                    tempi = wr * output[j + 1] + wi * output[j];
                    output[j] = output[i] - tempr;
                    output[j + 1] = output[i + 1] - tempi;
                    output[i] += tempr;
                    output[i + 1] += tempi;
                }
                wr = ((wtemp = wr) * wpr - wi * wpi) + wr;
                wi = wi * wpr + wtemp * wpi + wi;
            }
            mmax = istep;
        }

        for (i = 0; i < n; i++) {
            output[i] = output[i + 1];
        }
        return output;
    }
    /*public static double[] fftReal(double input[], int n, int isign) throws IOException{
    int i, i1, i2, i3, i4, n2p3;
    double theta,wtemp,wpr,wpi,wr,wi,output[];
    float c1 = 0.5f, c2, h1r, h1i, h2r, h2i;
    output = new double[n + 1];
    theta = Math.PI/n;
    if(isign == 1){
    c2 = -0.5f;
    output = SpectrumProcess.fft(input,n,1);
    }else{
    for(i = 0; i < 2*n; i++){
    output[i] = input[i];
    }
    c2 = 0.5f;
    theta = -theta;
    }
    wtemp = Math.sin(0.5*theta);
    wpr = -2.0*wtemp*wtemp;
    wpi = Math.sin(theta);
    wr = 1.0+wpr;
    wi = wpi;
    n2p3 = 2*n + 3;
    for(i = 2*n; i > 0; i--){
    output[i] = output[i - 1];
    }
    for(i=2; i<=n/2; i++){
    i4 = 1+(i3=n2p3-(i2=1+(i1=i+i-1)));
    h1r=(float)(c1*(output[i1]+output[i3]));
    h1i=(float)(c1*(output[i2]-output[i4]));
    h2r=(float)(-c2*(output[i2]+output[i4]));
    h2i=(float)(c2*(output[i1]-output[i3]));
    output[i1]=h1r+wr*h2r-wi*h2i;
    output[i2]=h1i+wr*h2i-wi*h2r;
    output[i3]=h1r-wr*h2r+wi*h2i;
    output[i4]=-h1i+wr*h2i+wi*h2r;
    wr=(wtemp=wr)*wpr-wi*wpi+wr;
    wi=wi*wpr+wtemp*wpi+wi;
    }
    for(i = 0; i < 2*n; i++){
    output[i] = output[i + 1];
    }
    if(isign == 1){
    output[0] = (h1r=(float)output[0])+output[1];
    output[1] = h1r - output[1];
    }else{
    output[0] = c1*((h1r=(float)output[0])+output[1]);
    output[1] = c1*(h1r - output[1]);
    output = SpectrumProcess.fft(output,n,-1);
    }
    return output;
    }*/
}
